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SUMMARY 


This  report  describes  a computational  method  called  "sharp  discontinuity 
tracking"  which  is  a new  approach  for  accurately  predicting  the  flow  field 
arising  in  explosion  problems.  This  work  was  performed  by  the  author  who  is  a 
member  of  the  Applied  Mathematics  Branch  of  NSWC.  The  guidance  and  advice  of 
Dr.  Jay  Solomon  of  NSWC  was  invaluable  and  is  gratefully  acknowledged. 
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1.  INTRODUCTION 


Many  problems  of  Interest  Involve  the  propagation  of  blast  waves  arising  from 
explosions,  and  the  interaction  of  these  waves  with  boundaries.  In  principle,  a 
computational  solution  of  this  fluid  dynamics  problem  should  provide  a cost- 
effective  alternative  to  expensive  physical  experimentation.  In  reality,  the 
computational  problem  is  an  extremely  difficult  one  which  involves  the  propagation 
and  interaction  of  various  discontinuities,  both  shock  waves  and  contact  disconti- 
nuities (e.g.,  interfaces  between  different  materials),  in  more  than  one  space 
dimension.  It  is  of  practical  interest  to  get  accurate  results  for  the  shock 
pressure  profile  (especially  the  peak  pressure)  and  also  to  correctly  predict  the 
motion  of  material  interfaces  which  significantly  affect  the  flow  field. 

Finite  difference  methods  currently  in  use  (both  Lagrangian  and  Eulerian)  for 
these  explosion  problems  use  either  explicit  or  implicit  artificial  viscosity  to 
spread  or  smear  shocks  over  several  computational  mesh  points.  This  falsifies  the 
peak  shock  pressure  and  impulse.  To  obtain  sufficient  accuracy,  it  is  necessary 
to  use  a very  large  number  of  mesh  points,  which  becomes  prohibitive  for  complicated 
multidimensional  problems.  While  pure  Lagrangian  methods  can  follow  material 
interfaces  well  and  can  allow  areas  of  fine  resolution  to  move  with  the  fluid,  they 
ultimately  tend  to  break  down  because  severe  material  deformation  leads  to  distorted 
computational  meshes.  Pure  Eulerian  methods  can  handle  distortion  and  internal  slip 
but  produce  diffusion  of  material  interfaces  and  do  not  readily  allow  localized 
resolution. 

This  research  is  aimed  at  devtloping  a new  numerical  method  called  "sharp 
discontinuity  tracking"  which  will  follow  shocks  and  contact  discontinuities 
accurately  without  introducing  smearing  or  computational  oscillations.  While  the 
method  is  specifically  designed  to  handle  the  problem  of  determining  the  pressure 
history  under  water  due  to  a chemical  explosion  in  air  above  a water  surface 
(herein  called  the  "EAAW  problem"),  it  should  be  applicable  to  other  explosion 
dynamics  problems  of  interest.  This  method  allows  for  adequate  computational 
resolution  of  physically  small  or  shrinking  regions  which  contain  important 
physical  phenomena,  such  as  the  air  region  between  the  expanding  explosive  gases 
and  the  water  surface  in  the  EAAW  problem  just  described. 

The  present  method  is  based  partly  on  the  work  of  Solomon  et  al.  [References 
1,2]  who  developed  a computational  method  for  predicting  the  three-dimensional 
steady  supersonic  flow  field  over  reentry  vehicles  at  angles  of  attack.  Since  such 
a steady  flow  problem  is  hyperbolic  with  the  axial  direction  being  a time-like 
coordinate,  it  is  similar  to  the  two-dimensional  axlsymmetric  unsteady  flow  case 
considered  here. 


1.  Solomon,  J.  M. , et  al. , "A  Program  for  Computing  Steady  Inviscid  Three- 
Dimensional  Supersonic  Flow  on  Reentry  Vehicles,  Vol  I:  Analysis  and 
Programming,"  NSWC/WOL  TR  77-28,  11  Feb  1977. 

2.  Solomon,  J.  M.,  et  al.,  "Inviscid  Flowfield  Calculations  for  Reentry  Vehicles 
with  Control  Surfaces,"  AIAA  Journal.  Vol.  15,  No.  12,  Dec  1977,  p.  1742. 
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2.  OVERVIEW  OF  COMPUTATIONAL  METHOD 

A.  GENERAL  DESCRIPTION  OF  METHOD.  The  method  of  sharp  discontinuity  tracking 

Is  a finite  difference  method  for  obtaining  approximate  solutions  of  the  nonlinear 
system  of  hyperbolic  conservation  laws  (partial  differential  equations)  which 
describe  invlscld  fluid  dynamics.  The  proposed  methodology  is  schematically 
indicated  in  Figure  1.  It  is  distinct  from  the  usual  finite  difference  methods 
since  it  defines  explicit  surfaces  in  the  computational  field  which  represent  * 

physical  discontinuities,  both  shock  waves  and  contact  discontinuities.  A special 
system  of  finite  difference  equations  governs  the  motion  of  these  surfaces  and 
exactly  enforces  the  physical  boundary  conditions  which  must  hold  across  the 
discontinuities.  No  finite  differences  are  taken  across  these  surfaces.  Since  the 
discontinuities  are  perfectly  sharp,  interactions  can  be  treated  by  locally  exact 
methods. 

In  practice  it  proves  useful  to  perform  a time-dependent  coordinate  trans- 
formation to  obtain  a regular  mesh  in  which  the  contact  surfaces  are  always 
coordinate  lines  and  the  shocks  "float"  through  the  Eulerian  mesh.  The  regions 
of  relatively  smooth  flow  between  the  discontinuities  can  be  accurately  handled  by 
standard  second  order  accurate  finite  difference  methods. 

A much  more  detailed  description  of  the  method  is  given  in  Appendix  A. 

B.  TRACKING  OF  DISCONTINUITIES.  The  surfaces  of  discontinuity  represent 
physical  shock  waves  and  contact  discontinuities.  The  motion  of  these  surfaces  is 
determined  by  the  physical  conditions  which  must  hold  across  the  surfaces  and  is 
also  influenced  by  the  surrounding  fluids.  The  correct  conditions  linking  the 
values  of  flow  variables  on  either  side  of  the  surfaces  are  the  Ranklne-Hugoniot 
jump  conditions  for  a shock,  and  the  equality  of  pressure  and  normal  particle 
velocity  for  a contact  discontinuity.  The  surrounding  fluid  can  influence  the 
surface  only  in  certain  ways  determined  by  the  theory  of  characteristics  for  a 
system  of  hyperbolic  partial  differential  equations.  By  selecting  appropriate 
admissible  characteristic  compatability  relations  and  combining  them  with  the 
correct  physical  boundary  conditions,  one  obtains  a system  of  differential 
equations  which  governs  the  motion  of  the  surface  and  gives  the  flow  variables  on 
either  side  as  a function  of  time.  These  differential  equations  are  then  discre- 
tized to  give  a system  of  finite  difference  equations  for  advancing  the  discontinu- 
ity in  time. 

C.  INTERACTION  OF  DISCONTINUITIES.  The  surfaces  of  discontinuity  move 
relative  to  each  other  and  will  in  general  collide  and  interact.  Where  they  meet 
they  are  locally  plane  and  the  values  of  flow  variables  on  either  side  are 
explicitly  known.  The  interaction  then  reduces  to  an  algebraic  problem  of  finding 
a configuration  of  transmitted  and  reflected  discontinuities  which  satisfy  all  of 
the  necessary  conditions  in  the  infinitesimal  neighborhood  surrounding  the  point 
of  Intersection.  This  locally  exact  solution  is  explicitly  inserted  into  the 
computed  flow  field. 

D.  EXPECTED  ADVANTAGES  OF  METHOD. 

(i)  Greater  accuracy  for  individual  discontinuities  since  they  remain 
perfectly  sharp  (not  smeared).  Hence  more  accurate  peak  shock  pressures  and 
impulses. 
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FIGURE  1 OVERVIEW  OF  COMPUTATIONAL  METHODOLOGY 
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(11)  Fewer  mesh  points  needed  since  shocks  are  not  spread  over  many  mesh 
points.  Decreased  computing  time  and  less  storage  necessary.  Possibility  of  doing 
more  complicated  problems  than  otherwise  feasible. 

(iii)  More  accurate  treatment  of  interactions  since  they  are  handled  by 
locally  exact  methods. 

(iv)  Localized  resolution  possible  within  a basically  Eulerian  computation 
by  varying  coordinate  transformation. 

3.  RESULTS 

A.  ONE  DIMENSIONAL.  The  shock  tracking  method  has  been  used  successfully  at 
NSWC  [Reference  1]  and  elsewhere  [Reference  3]  in  treating  aerodynamic  problems, 
but  apparently  has  never  been  applied  to  multidimensional  explosion  problems.  Our 
first  efforts,  however,  were  to  develop  the  tracking  method  for  spherical  explosion 
problems  in  one  space  dimension  and  time.  This  type  of  test  problem  is  ideally 
suited  for  the  standard  Lagrangian  finite  difference  approach  (as  contrasted  with 
our  methods  which  are  essentially  Eulerian)  since  mesh  distortion  is  no  problem  in 
one  dimension.  Furthermore,  there  are  special  features  of  the  one  dimensional 
problem  which  could  be  exploited  to  obtain  increased  computational  accuracy. 

However,  since  our  goal  is  to  develop  methods  for  multiple  space  dimensions,  we 
did  not  employ  any  techniques  in  the  one  dimensional  testing  which  wouldn't  be 
readily  extendable  to  more  space  dimensions.  The  following  cases  have  been 
computed: 

(i)  Spherical  Explosion  of  Pentolite  in  Water.  The  sequence  of  physical 
events  is  described  in  detail  in  [Reference  4].  Computational  results  for  peak 
shock  pressure  vs.  radial  distance  obtained  by  the  present  method  are  compared 
(Figure  2)  with  those  obtained  by  Sternberg  and  Hurwitz  [Reference  5].  These  latter 
results  are  known  to  agree  well  with  experiment.  It  is  noted  that  Sternberg  used 
a sharp  shock  algorithm,  but  his  method  is  not  readily  extendable  to  multiple 
dimensions.  He  also  took  advantage  of  the  fact  that,  in  one  space  dimension,  the 
initial  rarefaction  which  moves  back  into  the  explosive  gases  can  be  explicitly 
inserted  into  the  computation.  As  this  cannot  be  done  in  more  than  one  dimension, 
this  possibility  was  not  exploited  by  the  present  method.  Nevertheless,  agreement 
between  the  two  computational  results  is  excellent  (points  of  comparison  other  than 
those  shown  are  equally  close) . 


3.  Moretti,  G. , "Floating  Shock  Fitting  Technique  for  Imbedded  Shocks  in  Unsteady 
Multidimensional  Flows,"  in  Proceedings  of  the  1974  Heat  Transfer  and  Fluid 
Mechanics  Institute,  1974,  p.  184. 

4.  Sternberg,  H.  M.  and  Walker,  W.  A.,  "Calculated  Flow  and  Energy  Distribution 
Following  Underwater  Detonation  of  a Pentolite  Sphere,"  Physics  of  Fluids, 

Vol.  14,  No.  9,  Sep  1971,  p.  1869. 

5.  Sternberg,  H.  M.  and  Hurwitz,  H. , "Calculated  Spherical  Shock  Waves  Produced 
By  Condensed  Explosives  in  Air  and  Water,"  Proceedings  of  the  Sixth  Symposium 
(International)  on  Detonation,  Coronado,  Calif.,  Aug  1976,  p.  528. 


8 


162 


SHOCK  LOCATION  IN  CHARGE  RADII 

FIGURE  2 COMPUTED  SHOCK  PRESSURE  VS.  SHOCK  LOCATION, PENTOLITE  SPHERE  CENTRALLY 

DETONATED  IN  WATER  (X  INDICATES  STERNBERG/HURWITZ  RESULT  - OTHER  VALUES 
EQUALLY  CLOSE) 
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(il)  Spherical  Explosion  of  Pentolite  in  y = 1.4  Air.  This  problem  is 
more  difficult  than  the  explosion  in  water  case  because  the  initial  rarefaction 
moving  back  into  the  explosive  gases  is  considerably  stronger.  Shown  in  Figure  3 
is  a comparison  of  the  computational  results  for  peak  shock  pressure  vs.  radius 
obtained  with  the  present  method  (40  mesh  points  in  explosive  gas,  15  in  air  when 
main  air  shock  at  20  charge  radii),  those  obtained  by  Sternberg  and  Hurwitz 
[Reference  5],  and  a ccrve  fit  to  experimental  data  derived  by  Goodman  [Reference  6]. 
Agreement  is  generally  good,  despite  the  fact  that  the  present  method  does  not  » 

explicitly  insert  the  known  one  dimensional  rarefaction  solution  (as  does  Sternberg). 

The  pressure  profiles  as  a function  of  radius  at  selected  times  is  shown  in 
Figure  4,  where  the  qualitative  nature  of  the  interface  movement  and  secondary 
shock  propagation  is  seen  to  be  correct. 

B.  TWO  DIMENSIONAL  AXISYMMETRIC.  In  two  space  dimensions,  the  development  of 
the  method  is  incomplete,  so  that  the  entire  EAAW  problem  cannot  yet  be  handled. 

However,  some  preliminary  results  for  a spherical  explosion  in  air  above  a water 
surface  have  been  obtained.  These  results  are  shown  in  Figures  5 and  6 for  a 
centrally  detonated  explosion  of  a pentolite  sphere  located  30  charge  radii  above 
the  water  surface.  The  computation  used  an  extremely  coarse  finite  difference 
mesh  composed  of  about  500  mesh  points  (approximately  25  points  radially  on  each 
of  20  6 = constant  lines;  see  Appendix  A) . The  2-D  axisymmetric  computation  was 
initialized  with  the  1-D  computational  results  at  the  instant  before  the  spherical 
shock  hits  the  water  surface.  Figure  5 shows  the  main  shock  position  at  this 
instant  and  at  some  later  times.  Figure  6 shows  the  pressure  vs.  radius  at  these 
same  times  along  the  line  0=0  (straight  down  below  the  charge). 

No  comparison  with  experiment  is  yet  possible  because  the  computation  has  not 
been  carried  out  sufficiently  lar  in  time.  Further  results  cannot  be  obtained 
until  an  analysis  is  made  of  the  various  types  of  interactions  which  can  occur 
when  an  oblique  shock  hits  the  water  surface.  This  interaction  must  be  understood 
in  detail  because  the  present  computational  method  requires  an  explicit  local 
handling  of  discontinuity  interactions.  This  problem  is  more  difficult  than 
originally  anticipated  because  the  family  of  possible  interactions  appears  to  be 
very  complicated,  including  regular  refraction,  Mach  refraction,  and  refraction 
with  bound  and  fre  precursors  [Reference  7].  (So  far,  only  regular  refraction  has 
been  incorporated  in  the  present  method.)  Thus  this  seemingly  incidental  aspect 
of  the  overall  EAAW  configuration  is  in  fact  a difficult  fundamental  physical 
problem  in  and  of  itself.  Nevertheless,  the  preliminary  computational  results 
indicate  that,  once  these  local  interactions  are  understood,  the  present  computa- 
tional method  will  prove  successful  for  the  overall  explosion  dynamics  problem. 


6.  Goodman,  H.  J. , "Compiled  Free-Air  Blast  Data  on  Bare  Spherical  Pentolite," 
BRL  Rept.  1092,  Feb  1960. 

7.  Abd-El-Fattah,  A.  M.  , Henderson,  L.  F.,  and  Lozzi,  A.,  "Precursor  Shock  Waves 
at  a Slow-Fast  Gas  Interface,"  J.  Fluid  Mech. , Vol.  76,  1976,  p.  157. 
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FIGURE  4 PRESSURE  VS.  RADIAL  DISTANCE  AT 
VARIOUS  TIMES,  PENTOLITE  SPHERE 
CENTRALLY  DETONATED  IN  Y - 1.4  AIR 


FIGURE  5 MAIN  SHOCK  LOCATION  AT  VARIOUS  TIMES, 
PENTOLITE  SPHERE  CENTRALLY  DETONATED 
IN  AIR  30  CHARGE  RADII  ABOVE  FLAT 
WATER  SURFACE 


EXPLOSIVE  GAS H— AIR m-Um WATER 


FTCURE  6 PRESSURE  VS.  RADIAL  DISTANCE  AT 

VARIOUS  TIMES,  STRAIGHT  DOWN  BELOW 
CHARGE  (9=0),  PENTOLITE  SPHERE 
CENTRALLY  DETONATED  IN  AIR  30  CHARGE 
RADII  ABOVE  ELAT  WATER  SURFACE 
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APPENDIX  A 

DETAILED  DESCRIPTION  OF  THE  COMPUTATIONAL  METHOD 


This  appendix  contains  an  exposition  of  the  basic  mathematical  and  numerical 
techniques  which  underlie  the  sharp  discontinuity  tracking  method,  but  does  not 
give  a detailed  explanation  of  its  implementation. 


Al.  GOVERNING  EQUATIONS 


The  equations  of  inviscid  unsteady  flow  are  written  in  spherical  coordinates 
r,0,<f>  (FigureA-1)  for  an  axisymmetric  problem  (3/3<J>  + 0,  r i 0,  0 i 0 f n)  as  the 
following  system  of  weak  conservation  laws: 


where 


3U  3F 
3t  3r 


+ 


3G  . H 
30  + H 


0 


(A-l) 


F = r sin0 


■ pu 

2 

pu  +p 
puv 
u(pE+p)> 


-(2p+pv^)rsin0 
-r  (pcos0-puvsin0)/ 
0 / 


(A-2) 


Here  P is  the  density,  u is  the  r-component  of  velocity,  v is  the  0-component  of 
velocity,  e is  the  specific  internal  energy,  and  E = e + is  the  total 

specific  energy.  The  pressure  p is  given  by  an  equation  of  state  for  each  fluid 
(see  References  5,8  for  equations  of  state  for  pentolite,  air,  and  water)  in  the 
form 


p = p(p,e) 


(A- 3) 


8.  Walker,  W.  A.  and  Sternberg,  H.  M.,  "The  Chapman- Jouguet  Isentrope  and  the 
Underwater  Shock  Wave  Performance  of  Pentolite,"  Proceedings  of  the  4th 
Symposium  (International)  on  Detonation.  ACR-126,  U.S.  Gov't  Printing  Office, 
1965,  p.  27. 
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A "main"  shock  wave  is  defined  by  the  surface  r ■ s(9,t)  and  contact  dis- 
continuities are  given  as  r - cj(0,t),  i - 1,2, “'I.  For  the  explosion  in  air 
above  water,  1-2  and  c^  is  the  contact  surface  separating  the  explosive  gas 
products  from  air,  while  ci  is  the  air-water  surface  (plus  an  appropriate  closure 
for  large  6,  see  Figure  A-l).  The  main  shock  is,  in  general,  partially  in  air  and 
partially  in  water.  These  major  discontinuities  s and  c^  are  unknowns  which  must 
be  advanced  in  time  by  the  solution  procedure.  "Secondary"  shocks  and  contact 
discontinuities  are  allowed  to  be  smeared  out  by  the  finite  difference  scheme  used 
for  interior  mesh  points,  as  discussed  below. 

It  is  noted  that  the  choice  of  the  form  of  equations  system  (A-l)  eliminates  the 
need  to  specify  boundary  conditions  on  r * 0 and  on  the  symmetry  boundaries 
0 - 0,tt  since  the  unknown  U is  identically  zero  there.  This  is  consistent  with  the 
fact  that  by  symmetry,  u ■ v * 0 at  r * 0 and  v ■ 0 on  9 ■ 0,n.  However,  if  it  is 
desired  to  obtain  a solution  on  0 * 0,ir,  it  is  necessary  to  derive  another  set  of 
equations  which  hold  only  on  0 = 0 and  0 - it;  this  is  easily  done  by  dividing  (A-l) 
by  sin0  and  using  L'Hopital's  rule  to  evaluate  the  limits  0 -*•  0,w. 

A2.  COMPUTATIONAL  PROCEDURE 


The  computational  field  is  viewed  as  being  composed  of  regions  of  relatively 
smooth  flow  which  are  separated  from  each  other  by  explicit  discontinuities — either 
shocks  or  contact  discontinuities.  It  is  convenient  to  employ  a time  dependent 
transformation  which  always  maps  the  contact  discontinuities  c^,  which  are  single- 
valued, piecewise  smooth  curves  in  physical  space,  into  coordinate  lines  in 
computational  space  (Figure A-l).  It  will  be  seen  that,  for  the  explosion  in  air 
above  water  (EAAW) , this  has  the  advantage  of  preserving  computational  resolution 
in  the  air  region  c|  1 r 1 C2,  even  if  c 2 - cj  becomes  small. 

The  transformation  from  physical  to  computational  space  is  given  generally  as 

X - X(r,0,t) 

Y - Y(0)  (A-4) 

T - t 


and  the  governing  partial  differential  equations  (A-l)  are  expressed  in  the  trans- 
formed computational  space  (X,Y,T)  as 


where 


3U  3^ 

3T  3X  3Y 


ft 


0 


(A-5) 


ft  - XtU  + XrF  + X0C 


S-V 


(A-6) 


Once  a specific  choice  of  transformation  (A-4)  is  made,  the  quantities  Xt,Xr,Xg,YQ, 
and  their  partial  derivatives  are  given  explicitly  in  terms  o.f  s(9,t),  cj(0,t), 
and  their  partial  derivatives. 
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For  the  EAAW  problem,  convenient  choices  of  the  transformation  X are: 


gas  region: 

0 1 r t Cl(0,t) 


air  region: 
c1(6,t)  £ r < c2(6,t) 


water  region:  X - k..+k_[r-c9(0,t)] 

c2(0,t)  1 r 122 

so  that  r • ci(0,t)  is  mapped  into  X - 1,  and  r » C2(0,t)  is  mapped  into  X - kj. 
Then  for  a given  number  of  mesh  points  in  the  X-directlon  in  the  explosive  gas 
region,  the  choices  of  kj  and  k2  give  the  relative  computational  resolution  in 
the  air  and  water  regions  respectively,  A convenient  choice  of  the  transformation 
Y is  (in  inverted  form) 


where  w is  a parameter  which  gives  graded  resolution  in  physical  space  with 
finest  resolution  near  0 * 0 if  u > 1, 

lc 

A mesh  function  U . is  defined  on  the  uniform  computational  mesh 
ro.n.j 

for  X “ mAX  m * 0, 1, 2 , • • • ,M(j  ) 

j-1,2,3  Y - nAY  n - 0,1,2,***,N 

T - T +kAT  k - 0,1,2, •• • 
o 

where  m(l)AX  - 1,  m(2)AX  * kj  - 1,  NAY  « 1,  and  j*l,2,3  for  the  regions  of 
explosive  gas,  air,  and  water  respectively  (Figure  A-l) . The  main  shock  defined  by 
r ■ s(0,t)  in  physical  space  "floats"  through  the  uniform  mesh  and  is  represented 
by  special  mesh  points  which  are  located  on  Y * constant  lines. 

It  is  assumed  that  initial  data  (see,  e.g.,  Reference  8)  is  given  T * T0  for 
P,u,v,e,p,ci(0,to)  and  s(0,to)  on  the  spatial  mesh  defined  above.  In  order  to 
advance  the  solution  to  T0  + AT,  procedures  must  be  given  for 

(i)  "interior"  points  (1  < m 5 M(j)  -1,  1 * n i N - 1) 

(ii)  points  on  the  symmetry  boundaries  0 ■ 0,ir  (n  ■ 0,N) 

(iii)  shock  points 

(iv)  contact  discontinuity  points  (j  ■ 1,  m ■ M(l);  j * 2,  m ■ 1 and 
m * M(2) ; j ■*  3,  m • 1) 

(v)  interactions  between  a shock  and  a contact  discontinuity 
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The  solution  at  all  points  Is  advanced  in  time  using  finite  difference  predictor- 
corrector  methods.  For  the  interior  points  (i),  the  equations  (A-5)  are  discretized 
and  standard  MacCormack  [Reference  9]  schemes  (in  this  case,  predictor  forward 
differencing  and  corrector  backward  differencing  in  both  X and  Y)  are  used  to 
advance  the  conservation  vector  U.  The  primary  flow  variables  are  easily 
recovered  from  the  definition  of  U and  the  equations  of  state.  Points  on  the 
symmetry  planes  (ii)  are  advanced  in  the  same  way,  only  the  governing  equations 
must  be  modified  (as  described  above)  before  discretization.  The  procedures  for 
advancing  the  remaining  points,  namely  those  at  the  shock,  the  contact  discon- 
tinuities, and  the  points  of  discontinuity  interaction,  are  described  later. 

The  MacCormack  finite  difference  scheme  was  chosen  for  its  second  order 
accuracy  and  its  good  "shock-capturing"  properties  (ability  to  reasonably  smear 
out  a shock  over  several  mesh  points),  but  any  other  predictor-corrector  method 
could  be  used.  For  very  strong  shocks,  such  as  the  secondary  shock  which  implodes 
on  the  origin  in  a spherical  explosion,  it  is  necessary  to  implement  a computational 
"filter"  (similar  to  that  in  Reference  10  ) to  limit  shock  steepening.  Any  other 
shock  capturing  "improvement"  may  also  be  included. 

Computational  stability  demands  that  AT  be  suitably  restricted.  A necessary 
(but  not  sufficient  [Reference  11] ) criterion  can  be  derived  using  a geometric 
argument  which  requires  that  the  domain  of  dependence  of  the  "linearized"  differen- 
tial equations  be  contained  in  the  domain  of  dependence  of  the  finite  difference 
equations.  When  this  criterion  is  used  with  a multiplicative  safety  factor  of  .9, 
no  instability  is  observed  in  actual  computation. 

A3.  ADVANCING  THE  MAIN  SHOCK  AND  CONTACT  DISCONTINUITIES 

A crucial  requirement  of  the  present  computational  method  is  to  provide  an 
appropriate  numerical  treatment  of  the  explicit  surfaces  which  represent  shocks  or 
contact  discontinuities.  The  known  physical  boundary  conditions,  namely  the 
Rankine-Hugoniot  jump  conditions  and  the  equality  of  tangential  velocity  across  a 
shock,  and  the  equality  of  pressure  and  normal  velocity  across  a contact  dis- 
continuity, must  hold  exactly.  However,  these  boundary  conditions  alone  are  not 
sufficient  to  provide  equations  for  advancing  the  surface  in  time,  but  must  be 
considered  along  with  appropriate  information  obtained  from  the  governing  PDE 
system  (A-l).  From  the  theory  of  characteristics  for  hyperbolic  systems  of  PDE's, 
one  can  derive  characteristic  comparability  conditions  (see  Reference  12)  which 
are  associated  with  certain  bicharacteristic  directions.  Admissable  characteristic 


9.  MacCormack,  R.  W. , "The  Effect  of  Viscosity  in  Hypervelocity  Impact  Cratering," 
AIAA  Paper  No.  69-354,  1969. 

10.  Harten,  A.,  "The  Artificial  Compression  Method  for  Computation  of  Shocks  and 
Contact  Discontinuities.  I.  Single  Conservation  Laws,"  Comm.  Pure.  Appl.  Math, 
Vol.  XXX,  1977,  p.  611. 

11.  Turkel,  E. , "Phase  Error  and  Stability  of  Second  Order  Methods  for  Hyperbolic 
Problems  I,"  Journal  of  Computational  Physics.  Vol.  15,  June  1974,  p.  226. 

12.  Courant,  R.  and  Hilbert,  D. , Methods  of  Mathematical  Physics.  Vol.  II, 
Interscience,  N.Y.,  1962. 
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relations,  that  ia,  those  associated  with  directions  pointing  away  from  the 
surface  in  the  negative  time  direction,  effectively  tell  the  surface  what  is 
happening  in  the  surrounding  fluid.  In  the  manner  introduced  by  Kentzer 
Reference  13)  and  used  by  Solomon  et  al.  (References  1,  2),  it  is  possible  to 
combine  the  physical  boundary  conditions  with  the  correct  characteristic  compara- 
bility relations  to  obtain  a system  of  PDE's  which  hold  only  on  the  surface.  A 
finite  difference  approximation  of  these  equations  then  provides  the  algorithm  for 
advancing  the  surface,  and  the  flow  properties  on  either  side,  in  time. 

(i)  CONTACT  D I SCONT INI' ITY . A contact  discontinuity  (c.d.)  r ■ c(6,t)  is  a 
material  surface  for  the  fluids  on  either  side  of  the  discontinuity.  Hence 


u + — v “ 0 
c 


(A-8) 


holds  at  a c.d.  (subscripts  t,6  indicate  partial  derivatives).  There  are  three 
admissable  characteristic  compatability  relations  which  hold  on  fach  si^e  (left 
and  right)  of  the  surface.  Two  relations  correspond  to  a streamline  and  one  to 
the  Mach  conoid.  These  compatability  relations  are  of  the  form 

3p  » l£.  ■ r 
It  " Ai  8T  Ri 


3T  ^°cqtang^ 


A2  3T  + A3 


rang)  " “ Jr  (r)  ' R2 

Dr + * It  (r)]  - *3 


(A-9) 


(A-10) 


(A-ll) 


It  is  emphasized  that  (A-9)-(A-ll)  hold  on  each  aide  of  the  surface,  and  hence 
represent  six  equations.  Rj.  and  R2  contain  derivatives  interior  to  the  surface 
(Y-derivatives)  and  undifferentiated  quantities  evaluated  at  the  surface,  while 
R3  contains  additionally  X-derivatives.  Aj_,  A2,  and  A3  contain  only  undifferenti- 
ated quantities,  and  D*  - 1 + • The  pressure  p and  the  normal  velocity 


q - — are  constant  across  the  c.d.,  while  the  density  p,  energy  e,  and 
norm  D 

c 

tangential  velocity  9tang  can  sustain  jumps. 

In  addition  to  the  characteristic  compatability  relations  we  have  the  equation 


3 _/fe\  JB  fft  ft  3c] 

3T  l c / c [3Y  c 3YJ 


(A-12) 


I3j  Kentzer,  C.  P. , "Discretization  of  Boundary  Conditions  on  Moving  Discontinuities, 
Lecture  Notes  in  Physics.  Vol.  8,  Springer-Verlag,  1971,  p.  108. 
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which  expresses,  through  the  chain  rule,  the  equality  of  cross  partials 


i£  . v le.  + k — 

9T  1 3T  2 3T 


(A-13) 


holds  on  each  side,  where 


K = 

1 \3p/ 


and  K„ 
e z 


are  given  by  the  appropriate 


equations  of  state.  A procedure  is  now  clear  for  advancing  the  c.d.,  since 
equations  (A-9)-(A-13)  provide  nine  equations  in  the  nine  unknowns  which  are  (the 
time  derivatives  of)  ct,ce,p,[p,e,qtang]Left,[p,e,qt;ang]Rlght. 

(ii)  SHOCK.  A shock  r = s(0,t)  is  a surface  across  which  the  Rankine-Hugoniot 
jump  conditions  hold,  namely 


constant 


(A-14) 


p + pQ  = constant 
v xnorm 


(A-15) 


e + p/p  + Qnorm/2  = constant 


(A-16) 


where  <lnonn  is  the  velocity  normal  to  the  shock  in  (r,0)  coordinates, 


D - 1 + 
s 


and  Q 


q - — is  the  relative  velocity  normal  to  the 
norm  D 

s 


moving  shock  surface.  Additionally,  the  tangential  velocity  is  unchanged  across 
a shock, 

q = constant  (A-17) 

tang 

Only  shocks  moving  into  undisturbed  fluid  are  considered  here.  For  this  case, 
there  is  only  one  admissible  characteristic  compatability  relation,  and  it 
corresponds  to  a Mach  conoid  on  the  disturbed  (high  pressure)  side  of  the  shock. 

The  compatability  condition  is  of  the  form 


, 3q 

lE.  + A 71°™ 
3T  4 3T 


(A-lf  3 


where  R4  contains  derivatives  in  the  X-direction,  derivatives  interior  to  the 
shock,  and  undifferentiated  quantities.  Differentiating  (A-14)-(A-16)  with  respect  to 

T and  solving  for  and  ■n°rm , it  is  possible  to  rewrite  (A-18)  as 


A — t _ A — 
5 3T  a6  3T 


(A-19) 
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where  A5  and  Ag  contain  only  undifferentiated  quantities.  The  equality  of  cross 
partlals,  expressed  through  the  chain  rule,  is  of  the  form 


(A-20) 


so  that  (A-19)aid  &-20)  provide  equations  for  advancing  sq  and  st.  It  is  then  easy 

to  get  the  advanced  value  of  q and  use  the  Hugoniot  conditions  (A-14)-(A-16)  to 

norm 

solve  for  all  properties  behind  the  shock. 


(iii)  COMMENTS . For  either  a shock  or  a c.d.,  we  have  a system  of  PDE's  for 
advancing  (in  time)  the  surface  and  the  flow  variables  on  either  side.  These 
equations  are  implemented  with  a predictor-corrector  finite  difference  method. 
Spatial  differencing  internal  to  the  surface  is  handled  in  the  usual  predictor- 
forward,  corrector-backward  way.  However,  differencing  in  the  X-direction  requires 
special  treatment  since  no  differences  can  be  taken  across  the  discontinuity. 

Thus  at  the  left  side  of  discontinuities  the  differencing  is  always  backward  (both 
predictor  and  corrector) , whereas  on  the  right  side  the  differencing  is  always 
forward.  This  procedure  is  only  first  order  accurate,  but  it  is  possible 
(Reference  14)  to  add  a correction  term  to  achieve  overall  second  order  accuracy. 


At  the  shock,  finite  differences  are  taken  on  a non-uniform  mesh,  but  again 
it  is  possible  to  secure  second  order  accuracy  by  adding  a correction  term. 
Additionally,  since  the  shock  is  moving  into  a previously  undisturbed  fluid  where 
there  are  no  mesh  points,  it  is  necessary  to  insert  new  mesh  points  behind  the 
advancing  shock.  This  is  done  by  quadratic  interpolation  of  the  conservation 
vector  U using  the  known  primary  flow  variables  at  the  shock  and  the  values  of  U 
at  the  mesh  points  along  Y » constant  lines  behind  the  shock. 


A4.  INTERACTIONS 


Since  the  sharp  discontinuity  tracking  method  treats  physically  important 
discontinuities  as  explicit  surfaces,  it  is  necessary  to  analyze  the  situation 
which  occurs  when  two  such  discontinuities  interact.  Considered  here  is  the  case 
of  an  air  shock  hitting  the  air-water  surface  obliquely.  It  appears  that  the 
qualitative  nature  of  the  resulting  physical  phenomena  depends  upon  the  incident 
shock  strength  and  the  angle  at  which  the  discontinuities  meet  [Reference  7],  For 
sufficiently  strong  shocks  and  sufficiently  small  incidence  angles  a,  a "regular 
refraction"  solution  exists  as  shown  in  Figure  A-2.  In  a coordinate  system  moving 
along  the  surface  with  point  A,  the  flow  is  pseudo-steady.  Assuming  locally 
constant  states  near  A and  locally  plane  discontinuities,  the  problem  is  reduced 
to  an  algebraic  one  of  finding  the  twelve  unknowns  (p,u,v,e,p  in  regions  3 and 
4,  B.y)  which  satisfy  twelve  nonlinear  equations  (Rankine-Hugoniot  jump  conditions 
and  constancy  of  tangential  velocity  across  each  of  R and  W,  state  equations  in 
regions  3 and  4,  contact  discontinuity  conditions  across  the  disturbed  water 
surface) . 


14.  Wf-rming,  R.  F.  and  Beam,  R.  M. , "Upwind  Second-Order  Difference  Schemes  and 
Applications  in  Unsteady  Aerodynamic  Flows,"  Proceedings  of  the  AIAA  2nd 
Computational  Fluid  Dynamics  Conference.  Hartford,  Conn.,  1975,  p.  17. 
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explicitly  tracked 
incident  shock 


W 

transmitted  shock 


FIGURE  A- 2 REGULAR  REFRACTION  OF  AIR  SHOCK  AT  AIR/WATER  INTERFACE 
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This  locally  exact  solution  gives  the  necessary  information  to  track  the 
transmitted  water  shock  W and  the  disturbed  air-water  interface.  At  present  the 
reflected  shock  is  allowed  to  be  smeared  out  by  the  MacCormack  finite  differencing, 
but  there  is  no  conceptual  difficulty  in  explicitly  tracking  R as  well. 
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